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Numerical simulations based on electronic structure calculations are finding ever growing ap- 
plications in many areas of physics. A major limiting factor is however the cubic scaling of the 
algorithms used. Building on previous work [F. R. Krajewski and M. Parrinello, Phys.Rev. B 71, 
233105 (2005)] we introduce a novel statistical method for evaluating the inter-atomic forces which 
scales linearly with system size and is applicable also to metals. The method is based on exact 
decomposition of the fermionic determinant and on a mapping onto a field theoretical expression. 
We solve exactly the problem of sampling the Boltzmann distribution with noisy forces. This novel 
approach can be used in such diverse fields as quantum chromodynamics, quantum Monte Carlo or 
colloidal physics. 
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Atomistic simulations in which the interactions are 
computed on the fly from electronic structure calcula- 
tion play an important role in modern science and have 
proven their relevance in many fields. However their com- 
putational cost is a severe limitation. In particular sim- 
ulating large systems has proven challenging due to the 
cubic dependence of the computation time on the num- 
ber of electrons. This has long since been recognized0,Q 
and a number of algorithms have been suggested that in 
principle lead to linear scaling HSHUHHSIIS M ° St 
are based on the possibility in semiconductors or insula- 
tors of localizing the electronic orbitals. Linear scaling is 
then achieved by neglecting interactions between faraway 
atoms. This approach however suffers from poor conver- 
gence and leads to errors that are not easy to control. In 
metals the wavefunctions cannot be localized and only a 
handful of methods have been suggested |3, llC| . All in 
all it can be stated that in spite of considerable progress 
performing linear scaling ab-initio simulations is still very 
challenging. 

Very recently we have proposed a new algorithm which 
scales linearly in all physical dimensions for semiconduc- 
tors and metals alike [ll|. Here we reformulate the algo- 
rithm as a field theory. Sampling the resulting action 
is done stochastically. We show that, in spite of the 
statistical noise present in the evaluation of the forces, 
exact sampling can be performed. Our way of solv- 
ing this problem is general and can also solve prob- 
lems of similar nature that are encountered in quan tum 
chr omo dynamic s Hii lplj | , quantum Monte Carlo|l4l| and 
colloidal phvsics|lfij. where also the interaction can only 
be determined stochastically. 

Let us start with the generic expression for the total 
energy in theories that can be formulated in an effective 



single particle form: 

JV 

E = 2j2d + Vr- (1) 

i=l 

The first term is the so-called band structure term given 
by the sum of the lowest N doubly occupied eigenvalues 
6i of a Hamiltonian H. For instance in density functional 
theory H is the Kohn and Sham Hamiltonian and V r cor- 
rects for double counting and accounts for the direct nu- 
clear nuclear interaction, while in tight binding and other 
semi-empirical approaches H is a Hamiltonian which de- 
pends parametrically on the atomic positions and V r is 
a pairwise repulsive energy. In either case V r can be cal- 
culated in O (N) operations, while the calculation of the 
band structure term has an apparent O (-/V 3 ) complexity 
and has been the limiting factor that has so far prevented 
simulating very large sy stems. 

Following ref. [lfjllTj we write the band structure term 
as the low-temperature limit of the grand canonical po- 
tential for independent fermions: 

= -|lndet (l + e^" 11 )) . (2) 

In Eq. J2J n is the electron chemical potential and it easy 
to see that lim^oo £1 = 2 Yli=i £ i ~ M-^ej where N e = 2N 
is the total number of electrons. We now factorize the 
operator in Eq. J2J as: 

P/2 

1 + e -^- H ) = [] (MjM,) (3) 

i=i 

where P is an even integer and M; = 1 + 
e «7r(2;-i)/P e ^(^-H) jjere we depart from ref. ^l| an d 
since is a positive definite operator we can follow 

the practice well known in lattice gauge simulations^^ 
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of writing its inverse determinant as an integral over a 
field <&i that has the dimension of the full Hilbert space 
in the formH^: 



det ( m!m ; 



/2?[*,]e" 



(4) 



Inserting the relation © into Eq. (J2J after having used 
Eq. Q we end up with the following expression for the 
grand canonical potential: 



4 P/2 /■ 



,'mjMi#i + const. (5) 



which is the promised field theoretical formulation. The 
quantities of physical interest like energy or force can all 
be calculated as derivatives of f2 relative to an appro- 
priate external parameter. For instance N e = — , 
and assuming that /3 _1 is so low that temperature effects 
on the electrons can be neglected £; band = 2 Yli=i e « = 
^ (/3f2) + fJiN e while the contribution to force on parti- 
cle I at position R/ coming from the band term is given 
by Fj and = - Vr,0. In taking the derivatives the con- 
stant in Eq. © vanishes and one is left with calculating 
expressions of the type: 
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Thus all relevant properties can be evaluated by sampling 
the P/2 distributions e -*i M i M '*'. 

So far no approximation has been made and no com- 
putational advantage has been gained either. In order 
to make further progress we must take advantage of the 
fact that in M; the operator e"p^ _H ^ appears and that 
P can be taken to be sufficiently large for suitable ap- 
proximations to the exponential operator to be accurate. 
In ref. [jjj we used as basis set a grid in real space and 
a Trotter decomposition was the natural approximation 
to use. Here we will apply our method to a tight binding 
Hamiltonian and simply use a high-temperature expan- 



M, = l + e' r < 2, - 1 >/ p 



| 0* - H) 




(7) 



In such a manner the operator M; has the same sparsity 
of H, a fact which will eventually lead to linear scaling. 
We have also considered the possibility of using higher 
order expressions for the exponential operator in Eq. Q 
0. This leads to a smaller P value at the cost of a 
less sparse Mi- Since eventually we want to exploit par- 
allelism as much as possible we have preferred to use 
the expression (J7J). Note that no assumption has been 
made on the energy spectrum or on the local character 
of the wavefunctions and therefore our method will be 
valid both for metals and non-metals. It simplifies the 
calculation of the properties of the system if in Eq. © 
we use the expression 



— 2P« M < 
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which has an accuracy compatible with Eq. 0. In 
Eq. © O a can be /?, (fj, - H), or -/5V Rj H for A = fM, 
(3 or R/ respectively. A standard approach to sampling 
g-fjMjM,*! j g £ draw a sequence of normal distributed 
random numbers St/ and compute 3>z solving the equa- 
tions M/4>/ = SP/. Since M/ is sparse, this equation can 
be solved in O (N) operations using for instance a bi- 
conjugated gradient method, 19.]. It can be easily shown 
that this approach is equivalent to the stochastic inver- 
sion method advocated in ref. 111. Here we have decided 



instead to sample e 
ics. 



using a Langevin dynam- 



(9) 



where the components a of the white random noise vector 
obey the relation 



(efW(t)> = 2m ne 8(t) 



(10) 



(8) 



In this way we circumvent the need to invert the ma- 
trix M;. This solves the problem that for a tight bind- 
ing Hamiltonian as opposed to the Hamiltonian used in 
ref. 0] we were not able to find good preconditioners. 
Furthermore since different M/'s have different eigen- 
value spectra one can choose the to/ so as to achieve the 
optimal sampling speed in each I channel. This problem 
is particularly serious for metals where M i=; p can have 
eigenvalues close to zero, which leads to a much slower 
sampling speed [llj. Finally since we will use a Langevin 
sampling also for the ions it is pleasing to use the same 
sampling methodology for electronic and ionic degrees of 
freedom. 

Inevitably the interatomic forces that are calculated by 
this procedure will be affected by a statistical error. This 
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will prevent us from using these forces to perform en- 
ergy conserving MD calculations. However we will show 
that sampling of the Boltzman distribution is still pos- 
sible. Similarly to what was done for the electronic de- 
grees of freedom we sample the ionic configurations with 
a Langevin equation: 

MR/ = Fj - 7/R/ + Sj (11) 

where the random force obeys the relations 

(Bj(O)Bj(f)) = 6k B TM 7l 5(t) (12) 

and 

(F/(0)B/(i)} = . (13) 

From the electrons' Langevin dynamics we do not get 
the exact forces Fi but an approximation Ff = Fi + 
that is affected by a statistical error and therefore 
there is in principle no guarantee that correct Boltzman 
averages are obtained from the solution of Eq. (|llfl . But 
let us assume, and we will show later in a realistic case 
that this assumption is indeed justified, that also af is 
a white noise obeying 

(Ef(O)Ef(t)) 6k B TM-ffS(t) and 
(Fj(O)Bj(t)) - . (14) 

In this case the noise Bf simply adds to 3i and if we 
modify Eq. (|llfl so as to read 



MR/ =F/ - (7/ 



■7f) 



R/ + 3/ + Sf 



(15) 



we recover a Langevin equation whose trajectories can 
still be used to obtain a Boltzman sampling. 

At first sight it would appear that we are defeating our 
object since in general we know F/ + but not each 
term individually. However we can determine 7/ by vary- 
ing it until the equipartition theorem ^~MR|^ = ^k&T 

is satisfied. With this choice the sampling will be cor- 
rect and noisy forces can be used in the sampling. It is 
important to note that after the correct value for 7^ is 
determined it has to be kept constant during the simula- 
tion. With this procedure one can exactly calculate static 
observables within the framework of Langevin dynamics 
without knowing the exact force. This is at variance with 
ref. [2(lj where explicit assumptions on the noise have to 
be made. 

It remains now to show that the assumptions made 
above are correct and that the resulting scheme is ac- 
curate. As a test system we have chosen silicon in the 
diamond and in the liquid phase. We use the empirical 
tight binding interaction for silicon described in |2l| . We 
checked that a value of P = 200 is sufficient for the ap- 
proximation of Eq. Q to be valid. The Langevin dynam- 
ics parameters used were 5t e = 1, 7 e <5i e = 1/20 where 5t e 
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FIG. 1: a) Distribution of the kinetic energy for a sim- 
ulation of liquid silicon with 64 atoms at 3000K (cir- 
cles). Maxwell distribution (line). b) Autocorrelation 

function of the noise (Hf(O)Ef(i)) / ^ (F^ nd (0)) 2 ^ (line). 
Cross-correlation function of the noise and the exact force 
(F/(0)Hf'(t))/^(F5 >and (0)) 2 ^ (crosses). The correlation 

functions are normalized to the exact ionic force stemming 
from the band energy. All results are calculated for a system 
64 atoms of liquid silicon at 3000K. 



is the discretized inte gra tion time step and the algorithm 
of Ricci and Ciccotti |22j has been used throughout. The 
masses mi are adjusted such that the average force fluc- 
tuations are (m|M ; * ; ) & < 0.025. This adjustment 

of the masses is a crucial ingredient to obtain a good 
performance since if the masses were set to be equal for 
all values for I the fastest time scales of the field $p/2 
would be more than one order of magnitude larger than 
that of $1. The integration time step for the ion dy- 
namics is St = 1 fs. After each ionic displacement we let 
the <!>/ evolve under the action of the new M^M; until 
the distribution is equilibrated. The time needed for the 
equilibration is problem dependent and can be measured 
by looking at the correlation function (<&;(0)<&;(i;)). In 
the present case we make the rather conservative choice 
of running the electronic Langevin equation for 100 time 
steps. The *i thus obtained are used to calculate the 
ionic forces for the next integration step. The chemical 
potential is continuously adjusted such that the number 
of electrons fluctuates around the desired value. 

We first consider the case of 64 Si atoms in a periodi- 
cally repeated cell of length 10.86 A at a temperature of 
3000 K where the system is metallic. Using the proce- 
dure described above we find that if we take 7/ = ^ fs" 1 
we need to add a correction due to the noise in the forces 
j f = fs" 1 which fixes the average kinetic energy of 
the particles at the desired value. In Fig. ^ a) it is also 
seen that not only is the average energy correct but also 
that its fluctuations follow Maxwell distribution. In this 
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FIG. 2: (color) Pair correlation functions for liquid sili- 
con (3000K) (a) and crystalline silicon in the diamond phase 
(300K) (b). The red dots show the results from our new 
method and the lines are calculated using standard diagonal- 
isation of the tight binding Hamiltonian. 



algorithm is trivially parallelizable and therefore we ex- 
pect the situation to be much more favorable in terms of 
wall clock time on massive parallel platforms. Our code 
scales essentially linearly with P while the diagonalizion 
codes, although highly developed, are not amenable to 
similarly efficient parallelization. For instance Shellman 
et al.[23j have estimated that parallelization of a tight 
binding simulation on 1000 atoms loses efficiency when 
distributed over more than 8 processors. In addition we 
can gain even more if we organize our calculation hier- 
archically and distribute the most expensive part of the 
calculation, the matrix multiplication M|M/$/, to sev- 
eral processors. 

Finally we stress once again that the interest in our 
work transcends the field of linear scaling algorithms and 
offers an alternative to other MC methods with noisy 
estimators 0, |5(j . 

We thank A. Laio for his critical reading of the 
manuscript. 



small system it is possible to calculate the correct forces 
on the ions and check the statistical properties of H f . It 
is seen from Fig. Q] b) that Eq. (|14fl is satisfied to a very 
good approximation and that the correlation between 
the exact forces and the noise is minimal. Furthermore 
(Hj (0)Hj (£)) is strongly localized in time and its behav- 
ior is well fitted by two exponentials whose decay is much 
faster than the fastest ionic motion time scale. Therefore 
its net effect is similar to that of a delta function and 
in fact if we approximate (Hj (0)Hj (t)) as a delta func- 
tion whose strength is given by its integral we find an 
estimate for 7^ « fs _1 in good agreement with the 
empirical determination. Similar analysis conducted on 
semiconducting Si leads to a similar conclusion, namely 
that the effect of the noise in the calculation of the forces 
can be simply accounted for by adding the damping co- 
efficient 7j. In the solid phase we find jj = 555 fs -1 
which is similar to the result for the liquid. Therefore we 
expect that, with some adjustment phase transitions can 
be studied with our method, at least in this case. 

In Fig. we compare the pair correlation functions 
g (r) calculated with noisy forces and those evaluated 
with a standard approach. It is seen that the agree- 
ment is excellent and that the use of noisy forces does 
not degrade the quality of the simulation. 

There are two issues that we would like to discuss 
briefly. One is the break-even point between standard 
calculations and the present linear scaling method. This 
is not easy to determine since it depends on the accu- 
racy required, which in our case is related to the number 
P used in the decomposition of the fermion determinant 
Eq. (PJ. For the case of Si we estimate on a single proces- 
sor the crossing point to be at about 500 atoms. In other 
cases larger values of P might be necessary, thus shift- 
ing the crossing point to larger systems. However our 
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